Downloaded from http://biorxiv.org/on September 18, 2014 



bioRviv 

f V beta 

THE PREPRINT SERVER FOR BIOLOGY 

A Scalable Formulation for Engineering Combination Therapies for 
Evolutionary Dynamics of Disease 

Vanessa Jonsson, Anders Rantzer and Richard M Murray 
bioRxiv first posted online November 7, 2013 

Access the most recent version at doi: http://dx.doi.org/10.1101/000075 



Creative The copyright holder for this preprint is the author/funder. It is made available under 
Commons a CC-BY-NC 4.0 International license. 
License 



Downloaded from http://biorxiv.org/on September 18, 2014 



A Scalable Formulation for Engineering Combination Therapies for 

Evolutionary Dynamics of Disease 

Vanessa Jonsson, Anders Rantzer and Richard M. Murray 



Abstract — It has been shown that optimal controller synthesis 
for positive systems can be formulated as a linear program. 
Leveraging these results, we propose a scalable iterative algo- 
rithm for the systematic design of sparse, small gain feedback 
strategies that stabilize the evolutionary dynamics of a generic 
disease model. We achieve the desired feedback structure by 
augmenting the optimization problems with l\ and £2 regular- 
ization terms, and illustrate our method on an example inspired 
by an experimental study aimed at finding appropriate HIV 
neutralizing antibody therapy combinations in the presence of 
escape mutants. 

I. Introduction and Motivation 

A challenge inherent to the treatment of certain infectious 
and non-infectious diseases, such as HIV or cancer, is the 
risk that the pathogen or tumor will evolve away and become 
resistant to treatment methods that comprise the standard of 
care. Especially vulnerable to this phenomenon are treatment 
methods that involve exposing the disease population (such 
as viruses or cancer cells) to therapies targeting specific 
molecules involved in disease progression for an extended 
period of time. While these targeted therapies have the 
benefit of allowing physicians to tailor treatments to a 
patient's tumor cell population, they nonetheless establish an 
environment in which the occurrence of mildly drug resistant 
pathogens or tumor cells can develop an evolutionary advan- 
tage over those for which the therapy is targeted, leading to 
so called 'treatment-escape'. 

One of the solutions that has been proposed [1], [2] is 
the rational design of combination therapy, an approach that 
requires the identification of targets that are known to play 
a key role in disease progression. An example of a multi 
target therapy currently used for the treatment of chronic HIV 
infection is highly active antiretroviral therapy (HAART), 
which is comprised of a combination of antiretroviral drugs 
that target specific enzymes involved during different points 
of the infection cycle. 

Recent results by Rosenbloom, et al. [1] have been more 
quantitative in nature, modeling the evolutionary dynamics 
of HIV and showing through simulations how the effect of 
antiretroviral dynamics can determine HIV evolution and 
therapy outcome. The Michor lab [2] recently showed the 
effects of different erlotinib dosing strategies in the presence 
of pharmacokinetic fluctuations on the evolution of resistance 
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of non small cell lung cancer through simulations of a 
stochastic evolutionary dynamics model. 

Although these methods have provided some insight into 
the problem, the challenge of designing treatment protocols 
that prevent escape is one that has been more recently 
addressed by control theoretic methods. For cancer therapy, 
results in this spirit apply methods from optimal and receding 
horizon control [3], [4], as well as gain scheduling [5], to 
synthesize treatment protocols that are robust to parameter 
uncertainty, an inherent issue in all biological systems. In 
the context of HIV and antiretroviral therapy, the authors in 
[6] propose a discrete time formulation that allows for the 
design of switching therapy strategies to delay the emergence 
of highly resistant mutant viruses. 

In [7], we introduced a general algorithm that used an Woo 
approach for the principled design of targeted combination 
therapy concentrations that explicitly account for the evolu- 
tionary dynamics of a generic disease model. This algorithm 
was effective in generating robustly stabilizing controllers, 
however it suffered from an inherent lack of scalability symp- 
tomatic of semidefinite programming formulations. Here, 
we propose a scalable solution to the combination therapy 
problem by reformulating it as a second order cone program 
(SOCP). 

We also address the requirement that the synthesized 
controller be not only robust to unmodeled dynamics but 
also exhibit sparse structure and small feedback gains. This is 
motivated by the fact that the number of therapies commonly 
used in combination to treat a disease is often small while the 
number of potential usable therapies are often very large [8]. 
Targeted therapies such as small molecule drugs or antibodies 
exhibit a maximum effective concentration beyond which 
side effects are likely to worsen and no additional drug 
benefits are seen. 

In particular, through £ 1 and £ 2 regularization, we induce 
sparse structure in the feedback controller while bounding 
the magnitude of the feedback gains. This leads to a SOCP 
formulation of the combination therapy synthesis problem. 
The main contribution of this paper is a scalable algorithm 
for the systematic design of sparse, small gain feedback 
strategies to stabilize the evolutionary dynamics of a generic 
disease model. 

The article is structured as follows: In Section II, we re- 
call the extended quasispecies evolutionary dynamics model 
that encodes replication, mutation and neutralization and 
summarize relevant results in controller design of positive 
systems. In Section III, we present our C\ combination 
therapy synthesis algorithm. Section IV illustrates our al- 
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gorithm in the context of an HIV antibody therapy design 
problem previously studied in an experimental setting [9]. 
In this section we also compare performance and scalability 
properties of our C\ algorithm and the previously developed 
%oo algorithm [7] . Section V ends with concluding remarks 
and directions for future work. 



II. Preliminaries 



A. Notation 



R + denotes the set of nonnegative real numbers. The 
inequality X > 0,(X > 0) means that all elements of 
the matrix or vector X are positive (nonnegative). X y 0 
means that X is a symmetric and positive definite matrix. 
The matrix A £ M" xn is said to be Hurwitz if all eigenvalues 
have negative real part. Finally, the matrix is said to be 
Metzler if all off-diagonal elements are nonnegative. Define 
1„ to be the vector of all ones of dimension n. The induced 
matrix norm for a matrix M € 

||M|| p _ ind 



is 



sup 



\Mw\f 



w 



where \w\ p = (\ Wl \ p + + \w m \ p ) 1/p - Let G(s) = C(sl - 
A)^ 1 B + D be a rx m matrix transfer function. The 
induced norms of the corresponding impulse response g(t) = 
Ce At B + D5{t) are 



l|S||p-ind = SUp 



Iff * Hip 



for w € £PJ0, oo), given that g * w £ Cp[0, oo) is the 
convolution of g and w. Finally we refer to the oo-induced 
robust controller as the C\ controller as is customary in the 
robust control literature. 

B. Problem formulation 

The quasispecies model [10] was originally developed 
to describe the dynamics of populations of self replicat- 
ing macromolecules undergoing mutation and selection. We 
choose this model for its relative simplicity and its ability to 
capture the salient features of the evolutionary dynamics of 
a simplified generic disease model. In [7] we incorporated 
the effects of potential therapies into the basic quasispecies 
model, by defining a drug binding reaction, £ + x — ^> I ■ x 
— drug £ binds to self replicating macromolecule x with 
association rate Ka, giving a neutralized complex I ■ x. The 
extended quasispecies model for n mutants and m drugs, is 
written as: 

n m 

Xi = {nqu - d i )x l + ^ r iQikXk ~ ^2 ^ki^kXi (1) 



k^i 



k=l 



where X{ £ R + is the concentration of mutant i, 4 £ K + 
is the drug concentration (assumed to remain at constant 
concentrations throughout), and di are the replication 
and degradation rates, respectively, of mutant i, and is 
the probability that mutant k mutates to mutant i. Finally, 
i'ki = f(Kki) is a function of the association constant 
Kki for each neutralization reaction representing the rate at 



which a neutralizing macromolecule 4 neutralizes mutant 
i. The rates and ip^i can be viewed as replication and 
neutralization fitnesses of mutant i. When = 0, Vfe € 
the quasispecies dynamics are unstable. 
The following state space representation of equation (1) 
emphasizes the inherent feedback structure that arises from 
drug binding reactions: 



x = (A - *L)x 

z = Cx 



(2) 



with (i) A £ W lXn , with A i:j = r iqi3 > 0 V i ^ j and 
An — Tiqu — di, that encodes the mutation and replication 
dynamics; (ii) ^> £ M. nxnm a block diagonal matrix that 
describes the fitness of n mutants with respect to m different 
drugs, with diagonal elements ^ 



tyw) G Mi xm ; (iii) 



L = £ 



, with t = (4) £ M m , a block diagonal 



matrix that encodes the concentrations of drugs for all n 
mutants; (iv) C = [1„ L T ] T £ R( m " +1 ' xn ; and (v) w £ M'j 
an arbitrary positive disturbance. Note that £ E" xn is 
by construction a diagonal matrix. 

We set the regulated output z — [z\ Z2] T = Cx where 

z\ = xi -\ h x n minimizes the total virus population and 

Z2 = Lx serves as a proxy to minimizing the concentration 
of drugs needed to robustly stabilize the system. 

Remark 1: A is a Metzler matrix with off-diagonal entries 
that are several orders of magnitude smaller than the diagonal 
entries. This is due to the biological fact that mutation 
rates range from 10~ 5 — 10~ 9 mutations per base pair per 
replication cycle for reverse transcriptase to DNA replication. 

Letting G denote the closed loop system (2), the control 
task then becomes to design drug concentrations I by finding 
a controller L = (I eg I) that leads to a stable G satisfying 
HGHoo-ind < 7, for some robustness level 7 > 0. 

C. Linear programming controller synthesis for positive sys- 
tems 

Our previous work [7] provided a semi definite program 
(SDP) formulation to synthesize "Hoo controllers that stabi- 
lize the evolutionary dynamics of a generic disease model as 
described by the extended quasispecies model. One feature 
of our system was its internal positivity, which allowed us 
to restrict the storage function matrix used in the bounded 
real lemma to be strictly diagonal [11] [12]. Although this 
reduced the number of decision variables in the SDP, it 
was not enough to offset the inherent lack of scalability 
of semidefinite programming. Scalability is an important 
consideration when designing combination therapies where 
the number of possible mutants can be large, such as for the 
treatment of chronic HIV infection. Recent results [12], [13] 
on the synthesis of controllers for positive systems show that 
the design of structured static state feedback controllers for 
internally positive systems can be reformulated as a convex 
problem with methods that scale linearly with the number of 
non zero elements in the feedback matrix. In this section we 
provide a brief survey of the relevant definitions and results 
from [13]: 

Theorem 1: For the system: 
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x = (A + ELF)x + Bw 

z = (C + GLF)x + Dw { ' 

let V be the set of m x m diagonal matrices with entries in 
[0,1]. Suppose that A + ELF is Metzler and C + GLF > 0 
for all L £ T>. Let g(t) be the impulse response of 

G{s) = {C + GLF)[sI - (A+ ELF^B + D 

If the matrices B, D and F have non-negative coefficients, 
then the following two conditions are equivalent: 

1) There exists an L 6 T> with A + ELF Hurwitz and 

||.9||oo-ind < 7- 

2) There exists a £ e R™, fx € Rip with 

A£ + En + Bl< 0 
C£ + Gfi + Dl< 7 1 

If £, /i satisfy the linear constraints 2) then the stability and 
norm guarantees of 1) hold for every L such that [i = LF£. 

Input-output performance is characterized using induced 
norms which are determined by the closed loop system's 
static gain: 

Theorem 2: For anrxm transfer matrix G(s) — C(sl — 
A)- 1 B + D, let g(t) = Ce At B + DS{t) be the corresponding 
impulse response, where Ce At B > 0 for t > 0 and D > 0, 
with A Hurwitz. Then ||<?|| p _i n d = ||G(0)|| p _ind for p = 1, 
p = 2 and p = 00. 

The positive nature of the system allows us to restrict 
ourselves to linear storage functions, which in turn allows 
for sparse structure to be imposed on the feedback gain [12]. 
Our feedback gain L = I (g) £ is not only structurally con- 
strained to be block diagonal, but algebraically constrained 
as well, in that all block diagonal components must be equal. 
Unfortunately, there is no known convex reformulation for 
this additional constraint. 

D. Regularization for structured controller synthesis 

The biomedical justification for wanting a simple con- 
troller structure with small gains is twofold: first, the number 
of therapies that can be used simultaneously to treat a disease 
is often limited, and second to minimize side effects, it is 
desirable to keep the magnitude of drug concentrations small 
while being robust to evolution. These design specifications 
can be expressed with the use of regularization, a common 
technique used in machine learning and inverse problems for 
model identification [14], [15], [16], [17] and increasingly 
used for controller design [18], [19], [20], [21]. As such, 
we introduce l\ and £2 penalties in our design objective to 
promote controller sparsity and minimize controller gains. 

We combine these regularization techniques with con- 
troller synthesis results for positive systems and present an 
iterative algorithm that yields a suboptimal L\ controller. 
This formulation of the combination therapy problem allows 
the designer to explore explicit trade offs between closed 



loop performance, sparsity in controller structure and gain 
minimization. 

III. A Sub-optimal Li combination therapy 

CONTROLLER 

In this section, we address the aforementioned non- 
convexity of the optimal control problem by formulating an 
iterative algorithm for finding effective drug concentrations. 
Our main result addresses the issue of synthesizing a sta- 
bilizing controller subject to the constraints imposed by the 
quasispecies model (2), with acceptable robustness properties 
characterized in terms of its oo-induced closed loop norm. 

A. Initializing stabilizing controller 

We begin by synthesizing a stabilizing controller to use 
as an initial controller in our iterative algorithm, as noted 
in Remark 5. We recall a simple algorithm developed in [7] 
for the synthesis of a stabilizing controller, which admits a 
particularly simple formulation in light of the Metzler nature 
of A and the diagonal structure of ^L. 

Lemma 1: There exists e > 0 such that the solution to 
the convex program: 

minimize Nil I™ 

subject to (4) 

A d + el - -< 0 
L = I®£ 

is a stabilizing controller for system (2), where Ad £ R™ xra 
is a diagonal matrix such that (Ad)u = (A)u where A is the 
replication and mutation matrix in equation (2). 

Remark 2: The resulting LMI is strictly diagonal, and so 
reduces to an LP. 

Proof: Rewrite A = Ad + M where Ad is diagonal 
and M = {rriij} e R nxn , my = 0 for i = j and 
rriij > 0 for i 7^ j. By the Perron Frobenius theorem, 
there exists r > 0 such that the spectral radius p(M) — 
r < maxi 2~2 m n ■ Let e = maxj 2~2 rriij and rewrite M = 
el - {el - M). We note that -{el - M) -< 0. The closed 
loop dynamics are then given by A — ^L = Ad + el — {el — 
M)-^L -< Ad+el-^L -< 0, yielding the desired stability. 



B. A Suboptimal C\ combination therapy controller 

Observe that through a straightforward application of 
Theorem 2, with B = I, C = [1„ L T ] T , D = 0, E = 
—iff, F = I to system (1), solving the following non-convex 
program: 

minimize llCdloo + Aill^lli + A 2 ll^ll 2 
subject to 

Ax + Kx + 1 < 0 (5) 
K = *L 

L = I®£ 
x>Q 
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will yield a sparse combination of drug concentrations £, 
yielding an optimal oo-induced closed loop norm for appro- 
priately chosen regularizers Ai > 0, A 2 > 0 £ R. 

Remark 3: We can impose an additional constraint limit- 
ing the concentrations of candidate therapies. This is neces- 
sary with certain drugs that have maximum tolerated doses 
dictated by clinical trials. 

As mentioned earlier, there are no known convex refor- 
mulations of this problem due to the additional structure on 
L. As such, we suggest the following iterative algorithm, 
based on the convex programs (6) and (7), to find a stabi- 
lizing controller. For notation, let Y = Pz(x,j) denote an 
optimization problem P in which we optimize over x and 7 
leaving Z fixed and with solution Y. 
Program 1. Pl f (x,7) : 

minimize 7 
subject to 

Ax + ^>Lx + 1 < 0 (g) 

L = I®1, C=[l n L T ] T 

x>0 

7 > yc^iioo 

Program 2. P2 {X:XlM) (£,j) 

minimize 7 + Aill^lli + A2PII2 
subject to 

Ax + ^Lx + 1< 0 ( 7 ) 
L = I(g>e, C= [1„ L T ] T 
7 > ||C^||oo 
We now present our algorithm: 



Algorithm 1 Scalable Combination Therapy 

1) Set e > 0. 

2) Solve for initial stabilizing controller £': 

Solve (4) to obtain controller £ 0 - 

Set (x', 7 ) =Pl, 0 (x, 7 ). 
Set (£',7)=P2 (x , ;0 , 0) (£, 7 ). 

3) Find (A' x , A 2 , £) that minimize 7: 

V(Ai,A 2 ) e Ai x A 2 ,Ai,A2 £ R%, 
while 7' — 7 > e : 

Set (x',j) =Vl e (x,"f). 

Set (f )7 )=P2 (x , )Al>A2) (£,7). 

Set 7' = 7. 



Remark 4: The sequence of 7's defined by the iterative 
process in Algorithm 1 is non increasing by construction, 
and bounded below by 0, thus implying convergence. An 
initial stabilizing controller can always be found as shown 
in [7], and thus the algorithm can always be initialized. 
We therefore have that our algorithm always converges to 
a feasible, local minimum robustness value 7, generating a 
stabilizing controller for (2). 

Remark 5: In practice, we note that the C\ controller 
suffers from dependence on initial conditions and converges 



to local optima quickly, yielding a stabilizing controller with 
robustness properties that are not significantly different from 
the nominal controller. A solution to this is to iterate once 
through PI and P2, with Ai = A2 = 0 and initialize the 
algorithm with the resulting controller. 

Remark 6: Due to the presence of the £2 regularize^ P2 
and (5) are SOCPs and not linear programs. As it will be 
clearly demonstrated in the example in the next section, this 
is still more efficient than the SDP combination therapy 
algorithm from [7]. In addition, the second order cone 
constraint is only on the drug concentrations so it should 
have minimal effect on performance, given that the number 
of antibodies is small compared to the number of mutants. 

IV. HIV/Antibody Therapy Application 

Our results provide a principled approach to the design 
of antibody treatments for chronic infection with human 
immunodeficiency virus-1 (HIV-1). We illustrate this with an 
example motivated by experimental results of evolutionary 
dynamics of HIV-1 in the presence of antibody therapy 
obtained in [9]. 

A relatively recent discovery is that a minority of HIV- 
infected individuals can produce broadly neutralizing anti- 
bodies (bNAbs), that is, antibodies that inhibit infection by 
many strains of HIV [22]. These have been shown to inhibit 
infection by a broad range of viral isolates in vitro but also 
protect non-human primates against infection [22], [23], [24]. 
Recent experimental results conducted in the Nussenzweig 
lab at Rockefeller University have demonstrated that the use 
of single antibody treatments can exert selective pressure on 
the virus, but escape mutants due to a single point mutation 
can emerge within a short period of time [9]. Although 
antibody monotherapy did not prove effective, it was shown 
that equal, high concentrations of an antibody pentamix 
effectively control HIV infection and suppress viral load 
to levels below detection. The goal of this example is to 
demonstrate how our proposed algorithm offers a principled 
way to design combination antibody therapies that control 
HIV infection and prevent evolution of any set of known 
resistant mutants. In a realistic setting, the ability to do 
this relies on the knowledge of what resistant viruses may 
be selected for with single therapies, and so this algorithm 
would be most effective in conjunction with single antibody 
selection experiments. 

1 ) Model parameters: We consider a system of eighteen 
to thirty five HIV mutants with five potential antibodies to 
use in combination. Figure 5 lists the mutants that evolved 
from monotherapy experiments with their corresponding half 
maximal inhibitory antibody concentration (IC50) in fig /ml, 
as measured by the Nussenzweig lab in [9]. Antibodies 
3BC176, PG16, 45-46G54W, PGT128 and 10-1074 are po- 
tential combination therapy candidates. 

Although virus replication rates can vary considerably 
depending on the nature of the mutations a virus may 
undergo, we choose replication rates to be 0.5 (ml • day) -1 
for all mutants. We justify this selection by noting that 
escape mutants grew to be dominant mutants during selection 
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Full support controllers 

- ^''^-inrT 0 ' 5987 ^ =0, 1 723= ^ 
_L i . =0.67222/^0.17858 



2-sparse controllers 




| H controller 
■ L controller 



Fig. 1. The graphs depict the average of thirty simulations subject to 
random time invariant perturbations of 5.5 % in the plant dynamics found 
with the Haa and C\ combination therapy algorithms for evolutionary 
dynamics of the first 18 point mutants in Figure (5). (Left) Full support 
controllers synthesized with the pentamix of antibodies available: 3BC176, 
PG16, 45-46G54W, PGT128 and 10-1074 and (Right) Sparse controllers 
synthesized with only two antibodies 45-46G54W and PGT128. 



experiments and assume that replication rate variability due 
to mutations were negligible. 

The fitness function associated with the neutralization of a 
virus i with respect to an antibody j is a Hill function tpij = 
where n is the Hill coefficient, £j is the concentration 



of a given antibody j, and K, 



is the 



association constant for the virus/antibody binding reaction 

k 

lj + X; -^r £j ■ Xj, and k on and k 0 ff are the on and off 
reaction rate constants. Note that the association constant 
represents the fraction bound of antibody/virus complexes in 

3-IC50 

solution and that ify = 3r . + ;n( 2 )-ic50- ' * s f° un d by solving 
Equation (1) for one virus/antibody pair for the duration 
[£()>*/] = [0,3]. We simplify the Hill function by setting 
the Hill coefficient n = 1, as there is evidence that that 
antibodies do not bind cooperatively. Our algorithm yields 
antibody concentrations near zero and this yields the linear 
approximation tpij = jh-tij. In addition, the antibodies we 
consider in our example do not target the same epitope, in 
other words, do not bind competitively to the same sites on 
the virus, thereby reducing any coupling between antibody 
concentrations. 

2) Mutation process: The mutation rate for HIV reverse 
transcriptase is u — 3 x 10~ 5 mutations/nucleotide base 
pair/replication cycle, and the HIV replication cycle is ap- 
proximately 2.6 days. We approximate the rate of mutation 
for a particular amino acid mutation at a particular location 
to be — u(l — u) k = 1.443 x 10~ 6 per replication cycle, 
where k s» 3000 is the size of the genome in residues and 
n a = 19 is the number of amino acids that can be mutated 
to. Our model supports forward point mutations and two 
point mutations. We do not consider back mutations, as the 
probability of mutation is negligible. Units of concentration 
in number of viruses/ml or number of antibodies/ml are used 
for states, and time is measured in days. The standard volume 
is 1 ml. 

A. A comparison between T-Loo and L\ controllers 

We seek to provide a qualitative comparison between 
controllers synthesized using the scalable C\ and the Hoo 
algorithms. To do this, we adapt the formulation in [7] 




Fig. 2. The graph depicts the average (left) norm and oo — ind norm 
as a function of the sparsity of controllers found using either the if oo or 
L\ combination therapy algorithms. 



to include £ x and £ 2 regularization terms and solve the 
following non convex problem using our iterative algorithm. 



minimize 7 + Ai||£||i + A2PH2 
subject to 

A T cl X + XA ci + C T C X 



-< 0 



(8) 



X - 7 2 / 
A cl = (A- VL) 
C = [1 T L T } T 
L = I®£ 
X y 0, X diagonal 

We synthesized a nominal stabilizing controller using (4), 
a robust controller that minimizes the Hoo closed loop norm 
using (8), and a robust controller using (7) that minimizes 
the L\ closed loop norm for the evolutionary dynamics 
of the eighteen HIV point mutants listed in Figure 5. We 
found similar gains and robustness properties for both sparse 
and full controllers using either algorithm with the notable 
difference seen in computational time. Not surprisingly, the 
Ci algorithm has far superior performance, beating the 
runtime for the Woo synthesis algorithm by four orders of 
magnitude. 

We averaged thirty simulations of closed loop evolu- 
tionary dynamics subject to 5.5% random time invariant 
perturbations in the plant dynamics using both sparse and 
full support Woo and L\ controllers. The sparse controller 
found by the L\ algorithm performed better than the one 
found by "Hco algorithm, whereas the situation was reversed 
for the respective synthesized full support controllers. As 
previously mentioned, the motivation for generating sparse 
controllers for combination therapy is that number of ther- 
apies commonly used in combination to treat a disease is 
often limited for clinical reasons. Therefore, the potential 
for the Ci algorithm to synthesize controllers that are not 
only sparse but more robustly stable than Hoc algorithm is 
a desirable feature. 

Figure 2 shows the relationship between gain sparsity and 
both Hoo and oo-induced norms in the synthesis of Hoc 
and Ci controllers. Although closed loop Hoo norms remain 
constant with respect to sparsity for both controllers, the 
closed loop oo-induced norm decreases with sparsity with 
the suggesting that the C\ synthesis algorithm, as expected, 
finds better performing sparse controllers. 

We note that due to computational limitations, we were 
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not able to synthesize controllers using the Ti.^ algorithm 
for the full set of thirty five mutants from Figure 5. 

B. L\ controller synthesis, full HIV example 

We synthesized a nominal stabilizing controller using (4) 
comprised of an antibody pentamix (0.4687,0.7815, 0.6129, 
0.6279, 0.8831) ng/ml of (3BC176, PG16, 45-46G54W, 
PGT128, 10-1074), and a robust controller using (8) that 
consisted of antibody trimix (0.6891,0.6712,1.0706) /ml 
of (3BC176, 4546-G54W, PGT128). These both were gen- 
erated for the evolutionary dynamics of the full, thirty five 
HIV mutants listed in Figure 5. 

Both antibody pentamix (stabilizing) and trimix (robustly 
stabilizing) controllers have similar gains and based on a 
cursory first glance, one might be tempted to believe these 
have comparable robustness properties. Indeed, for some 
simulations of the closed loop dynamics subjected to 5% 
random time invariant perturbations in plant dynamics, the 
nominal controller is stabilizing as seen in Figure 3. This 
is qualitatively consistent with the experimental results done 
in by the Nussenzweig lab in [9]. It was shown that with 
weekly injections of equal concentrations of the antibody 
pentamix holding concentrations constant, had viral loads 
that remained below the limit of detection during an entire 
treatment course in mice. 

In [9] again, an antibody trimix of equal concentrations of 
3BC176, PG16 and 45-46G54W was suggested and experi- 
mentally shown to produce a decline in the initial viral load. 
However, a majority of mice in the experimental study had a 
viral rebound to pre-treatment levels, suggesting that in these 
cases, the virus had evolved mutations that were resistant to 
the trimix treatment. To compare the performance of our C\ 
synthesized controller with gains of (0.6891,0.6712,1.0706) 
fig/ml of (3BC176, 4546-G54W, PGT128) to the exper- 
imentally studied trimix, we chose equal concentrations 
of (3BC176, PG16, 45-46G54W), namely (1, 1, 1) ng/ml 
for the experimentally derived trimix. We found that even 
though total antibody concentrations were larger in our 
version of the experimental trimix, the robustly stabilizing 
controller synthesized by the L\ algorithm nonetheless per- 
formed overall better; the closed loop norms were HGjloo = 
0.2941 and ||G||co-ind= 0.6533 for the C\ controller versus 
IIGHoo^.26433 and || G ? || oo _i nd =0.74572 for the experimen- 
tal trimix. 

These simulations demonstrate that although many stabi- 
lizing solutions to the combination therapy problem exist, 
the best ones are found when design parameters such as a 
sparsity, limits on the magnitude of gains, and robustness 
guarantees are simultaneously considered. Experimentally 
searching for these combinations is infeasible as the number 
of potential therapies and possible concentrations to consider 
is experimentally intractable. We propose to guide these ex- 
perimental activities with our ability to design and synthesize 
combination therapy controllers. As such, one could generate 
a family of controllers based on "design specifications" 
tailored not only the (viral or cellular) composition of the 
disease, but to explore tradeoffs between number of therapies 
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Fig. 3. Sum of virus populations subject to random time invariant 
perturbations of 5% in the dynamics for 30 different simulations for 
(left) a stabilizing closed loop controller comprised of antibody pentamix 
(0.4687,0.7815, 0.6129, 0.6279, 0.8831) p,g/ml of (3BC176, PG16, 45- 
46G54W, PGT128, 10-1074) synthesized using the convex program (5) and 
(right) a robustly stabilizing closed loop controller comprised of antibody 
trimix (0.6891,0.6712,1.0706) pg/ml of (3BC176, 4546-G54W, PGT128) 
synthesized using the C\ combination therapy algorithm. 

used (sparsity), therapy concentrations (magnitude of the 
gain) and ability to support pharmacokinetic fluctuations 
(robustness to perturbations) and subsequently verify these 
experimentally. 

V. Conclusion and Future Work 

Leveraging recent results in positive systems, we pro- 
posed a scalable SOCP based iterative algorithm for the 
systematic design of sparse, small gain feedback strategies 
that stabilize the evolutionary dynamics of a generic disease 
model. Through the addition of l\ and £2 regularization 
terms to the objective function, we achieved the desired 
feedback structure. In future work, we plan to explore a 
principled integration of our methods with recent results on 
the robust C\ stability of positive systems [25]. In particular, 
we hope to explicitly account for model error introduced due 
to model linearization, parametric uncertainty and unmodeled 
dynamics due to drug interactions. 
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Fig. 4. Stabilizing gains found for nominal stabilizing controller, a robust controller using (8) that minimizes the 'Hoo closed loop norm and a robust 
controller using (6) that minimizes the C± norm of the closed loop system for evolutionary dynamics systems of the first eighteen HIV point mutants listed 
in Figure 5. 
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Fig. 5. IC50 values for the indicated antibodies on YU2 mutant viruses found in continuous antibody mono therapy experiments conducted by the 
Nussenzweig lab at Rockefeller University [9]. The trimix of antibodies is : 3BC176,PG16,45-46G54W, the penta-mix is 3BC176, PG16, 45-46G54W , 
PGT128 and 10-1074. Estimated two point mutations represent intermediary mutations needed for our model but not included in experimental results in 
[9]. The IC50 values were taken to be the maximum IC50 of both mutations. 
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